### This script will extract all static spatial data for all sites for BRT analysis
### Script was appended to add dem, slope, and aspect to 'final' dataset

# Load library SDMTools

library('SDMTools')

# Establish directories



# Read in XY data and microclimate data

xy = read.csv('/home1/99/jc152199/underlogdownscale/ul_xy.csv',header=T)

#Read in data frame with model data

micro = read.csv('/home1/99/jc152199/underlogdownscale/rawmicroclimatedata/ulagWITHOUTEmpiricalAirTemp.csv', header=T)

# Read in static ASCII files using a loop and assigning names from asciinames

asciis = list.files('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC', full.names=T)[-2]

asciinames = gsub('_WTplusbuffer_LatLong_WGS1984_250mres.asc','',basename(asciis))

i=1

for (ascii in asciis)

	{
	
	assign(asciinames[i],read.asc(ascii))
	
	cat('\n',asciinames[i],' Reading Completed\n',sep='')
	
	i=i+1
	
	}
	
### Done

# Extract data for static ASCII's

xy$d2c = extract.data(cbind(xy$long,xy$lat),coastdist)
xy$d2s = extract.data(cbind(xy$long,xy$lat),logplusonestreamdist)
xy$fpc = extract.data(cbind(xy$long,xy$lat),fpc)
xy$dem = extract.data(cbind(xy$long,xy$lat),dem)
xy$slope = extract.data(cbind(xy$long,xy$lat),slope)
xy$aspect = extract.data(cbind(xy$long,xy$lat),aspect)

# Merge xy with micro

micro_xy = merge(micro,xy,by=c('site'))

####################

### Needs to include data on log size

logsize = read.csv('/home1/99/jc152199/underlogdownscale/cjs_logsizes.csv',header=T)

#### Concatenate together year,month, and date to make a comparable integer

logsize$intstartdate = as.integer(paste(logsize$startyear,sprintf('%02i',logsize$startmon),sprintf('%02i',logsize$startday),sep=''))

### Do the same for micro_xy

micro_xy$intdate = as.integer(gsub('-','',micro_xy$date))

### Create some blank columns in micro_xy

micro_xy$logcirc = NA
micro_xy$loglen = NA
micro_xy$logden = NA
micro_out=NULL

#### Now match log data onto micro_xy

for (site in unique(logsize$point_ID))

	{
	
	### Subset to a single site

	sub_micro = micro_xy[which(micro_xy$site==site),]
	sub_log = logsize[which(logsize$point_ID==site),]
	
	### Basic match
	
	sub_micro$logcirc = sub_log$AvgOfCircumference[1]
	sub_micro$loglen = sub_log$Length[1]
	sub_micro$logden = sub_log$AvgOfPenetrometer[1]
	
	if(nrow(sub_log)>1)
	
		{
		
		sub_micro$logcirc[which(sub_micro$intdate>=sub_log$intstartdate[2])]= sub_log$AvgOfCircumference[2]
		sub_micro$loglen[which(sub_micro$intdate>=sub_log$intstartdate[2])] = sub_log$Length[2]
		sub_micro$logden[which(sub_micro$intdate>=sub_log$intstartdate[2])] = sub_log$AvgOfPenetrometer[2]
	
		}
		
	micro_out = rbind(micro_out,sub_micro)
	
	cat('\n',site,' - Completed\n')
		
	}
	
### Done with loop

### Calculate log volume

micro_out$logvol = (micro_out$logcirc/(2*pi))*(micro_out$logcirc/(2*pi))*pi*micro_out$loglen

#### Now need to intersect micro_out with BRT predictions of daily Tmax and daily Tmin

brt.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/brtpredsfinal/'

### Need to remove dates after 20090601 because there are no BRT estimates of temp for those days
#### Removes approximately 4000 records and completely removes 0 sites

micro_out = micro_out[which(as.numeric(gsub('-','',micro_out$date))<=20090601),]

#### Remove dates from February 29th (for which no BRT temperature estimates exist)

micro_out = micro_out[-which(as.numeric(substr(gsub('-','',micro_out$date),5,8))==0229),]

### Create some blank columns in micro_out and a blank object for binding within the loop

micro_out$BRTairmax = NA
micro_out$BRTairmin = NA

outdata = NULL

i=1

### Use a for loop to intersect with BRT temperature estimates

for (date in unique(gsub('-','',micro_out$date)))

	{
	
	### Subset to a single date
	
	tmo = micro_out[which(gsub('-','',micro_out$date)==date),]

	#### Identify BRT layers of interest
	
	maxbrt = read.asc.gz(paste(brt.dir,'maxgzip/',substr(date,1,4),'/tmax.',date,'.asc.gz',sep=''))
	minbrt = read.asc.gz(paste(brt.dir,'mingzip/',substr(date,1,4),'/tmin.',date,'.asc.gz',sep=''))
	
	### Intersect data
	
	tmo$BRTairmax = extract.data(cbind(tmo$long,tmo$lat),maxbrt)
	tmo$BRTairmin = extract.data(cbind(tmo$long,tmo$lat),minbrt)
	
	#### Bind data
	
	outdata = rbind(outdata,tmo)
	
	### Track position
	
	i=i+1
	
	#### Report progress
	
	cat('\n',(i/length(unique(micro_out$date)))*100,' - Percent Complete\n',sep='')
	
	}
	
### Close loop

### Write out outdata

write.csv(outdata, file=paste('/home1/99/jc152199/underlogdownscale/ModelDataforULBRT.csv',sep=''),row.names=F)
		
		
		
	









